/*
点到面的距离，在dir一侧为正，另一侧为负
点到直线距离
直线到直线距离
直线与平面交点
面与面的交线
两直线交点
三角形面积
p是否在三角形ps中
线段AB与三角形ps是否相交, p是交点
三角形与三角形是否相交
点p到线段AB距离
线段p到线段q不平行时的距离s, p q平行返回0不计算s
点投影到面上
三角形之间距离
三维凸包
3D 外接圆圆心
四面体外接球圆心
最小包围球
面切球

待测:
点在线段上
onSeg(){}
onSeg_Strict(){}


待加入:
将三维平面旋转至水平二维平面，使用二维计算，见convexHull2D
trans3D2D(){}
ftrans3D2D(){}
线段是否相交
int SinS(){}
点投影到线
getprojL(){}
点关于面镜像
reflect(){}

*/

#include <bits/stdc++.h>
#define rep(i,a,b) for(int i=(a);i<=(b);++i)
#define per(i,a,b) for(int i=(a);i>=(b);--i)
#define pb push_back
#define fi first
#define se second
using namespace std;

typedef long double db;
const db EPS = 1e-8;
const db INF = 1e100;
const db PI=acos(-1.0);
inline int sgn(db a) { return a < -EPS ? -1 : a > EPS; }
inline int cmp(db a, db b) { return sgn(a-b); }
struct P {
    db x, y;
    P(){}
    P(db _x, db _y) : x(_x), y(_y) {}
	bool operator<(P o) const { return cmp(x,o.x) != 0? cmp(x,o.x) == -1 : cmp(y,o.y) == -1; }
	bool operator>(P o) const { return cmp(x,o.x) != 0? cmp(x,o.x) == 1 : cmp(y,o.y) == 1; }
	bool operator==(P o) { return cmp(x,o.x) == 0 && cmp(y,o.y) == 0; }
    P operator+(P p) { return P(x + p.x, y + p.y); }
    P operator-(P p) { return P(x - p.x, y - p.y); }
    P operator-() { return {-x,-y}; }
    P operator*(db d) { return P(x * d, y * d); }
    P operator/(db d) { return P(x / d, y / d); }
    db dot(P p) { return x * p.x + y * p.y; }
    db det(P p) { return x * p.y - y * p.x; }
    db disTo(P p) { return (*this - p).abs(); }
    db alpha() { return atan2(y, x); }
    void read() { cin >> x >> y; } // 根据题目修改
	void output() { cout << x << "," <<y << '\n';}
    db abs() { return sqrt(abs2()); }
    db abs2() { return x * x + y * y; }
    P rot90() { return P(-y, x); } // 逆时针
	P rot(db a) { db c = cos(a), s = sin(a); return {c*x-s*y, s*x+c*y}; } // 绕原点 a 弧度
    P unit() { return *this/abs(); }
    int quad() const { return sgn(y) == 1 || (sgn(y) == 0 && sgn(x) >= 0); } // 判断是否在上半区
};

#define cross(p1,p2,p3) ((p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y))
#define crossOp(p1,p2,p3) sgn(cross(p1,p2,p3))

vector<P> convexHull(vector<P> ps) { // 相邻边夹角小于180
    int n = ps.size(); if(n <= 1) return ps;
    sort(ps.begin(), ps.end());
    vector<P> qs(n * 2); int k = 0;
    for (int i = 0; i < n; qs[k++] = ps[i++])
        while (k > 1 && crossOp(qs[k - 2], qs[k - 1], ps[i]) <= 0) --k;
    for (int i = n - 2, t = k; i >= 0; qs[k++] = ps[i--])
       	while (k > t && crossOp(qs[k - 2], qs[k - 1], ps[i]) <= 0) --k;
    qs.resize(k - 1);
    return qs;
}

struct P3{
	db x,y,z;
	P3(){}
	P3(db _x,db _y,db _z){x=_x,y=_y,z=_z;}
	P3 operator + (P3 k1){return P3(x+k1.x,y+k1.y,z+k1.z);}
	P3 operator - (P3 k1){return P3(x-k1.x,y-k1.y,z-k1.z);}
	P3 operator * (db k1){return P3(x*k1,y*k1,z*k1);}
	P3 operator / (db k1){return P3(x/k1,y/k1,z/k1);}
	db abs2(){return x*x+y*y+z*z;}
	db abs(){return sqrt(x*x+y*y+z*z);}
	db disTo(P3 o){return ((*this)-o).abs();}	
	P3 unit(){return (*this)/abs();}
	int operator < (const P3 k1) const{
		if(cmp(x,k1.x)!=0)return x<k1.x;
		if(cmp(y,k1.y)!=0)return y<k1.y;
		return cmp(z,k1.z)==-1;
	}
	int operator == (const P3 k1){ 
		return cmp(x,k1.x)==0&&cmp(y,k1.y)==0&&cmp(z,k1.z)==0; 
	}
	P3 det(P3 k2){ return P3(y*k2.z-z*k2.y, z*k2.x-x*k2.z, x*k2.y-y*k2.x); }
	db dot(P3 k2){ return x*k2.x+y*k2.y+z*k2.z; }
	void read(){ double k1,k2,k3; scanf("%lf%lf%lf",&k1,&k2,&k3);x=k1;y=k2;z=k3;}
	void output() { cout << x << "," <<y <<","<<z<< '\n';}
};

//P3 cross(P3 k1,P3 k2){return P3(k1.y*k2.z-k1.z*k2.y,k1.z*k2.x-k1.x*k2.z,k1.x*k2.y-k1.y*k2.x);}
//db dot(P3 k1,P3 k2){return k1.x*k2.x+k1.y*k2.y+k1.z*k2.z;}

//p=(3,4,5),l=(13,19,21),theta=85 ans=(2.83,4.62,1.77)
P3 turn3D(db k1,P3 l,P3 p){
    l=l.unit(); P3 ans; db c=cos(k1),s=sin(k1);
    ans.x=p.x*(l.x*l.x*(1-c)+c)+p.y*(l.x*l.y*(1-c)-l.z*s)+p.z*(l.x*l.z*(1-c)+l.y*s);
    ans.y=p.x*(l.x*l.y*(1-c)+l.z*s)+p.y*(l.y*l.y*(1-c)+c)+p.z*(l.y*l.z*(1-c)-l.x*s);
    ans.z=p.x*(l.x*l.z*(1-c)-l.y*s)+p.y*(l.y*l.z*(1-c)+l.x*s)+p.z*(l.x*l.x*(1-c)+c);
    return ans;
}
typedef vector<P3> VP;
typedef vector<VP> VVP;

//***********************?? 不懂
db Acos(db x){return acos(max(-(db)1,min(x,(db)1)));}

// 球面距离 , 圆心原点 , 半径 1
db Odist(P3 a,P3 b){db r=Acos(a.dot(b)); return r;}
db r; P3 rnd;
vector<db> solve(db a,db b,db c){
    db r=sqrt(a*a+b*b),th=atan2(b,a);
    if (cmp(c,-r)==-1) return {0};
    else if (cmp(r,c)<=0) return {1};
    else {
        db tr=PI-Acos(c/r); return {th+PI-tr,th+PI+tr};
    }
}
vector<db> jiao(P3 a,P3 b){
    // dot(rd+x*cos(t)+y*sin(t),b) >= cos(r)
    if (cmp(Odist(a,b),2*r)>0) return {0};
    P3 rd=a*cos(r),z=a.unit(),y=z.det(rnd).unit(),x=y.det(z).unit();
    vector<db> ret = solve(-(x.dot(b)*sin(r)),-(y.dot(b)*sin(r)),-(cos(r)-rd.dot(b)));
    return ret;
}
db norm(db x,db l=0,db r=2*PI){ // change x into [l,r)
    while (cmp(x,l)==-1) x+=(r-l); while (cmp(x,r)>=0) x-=(r-l);
    return x;
}
//************************
//  把一个面变成 ( 点 , 法向量 ) 的形式
pair<P3,P3> getF(VP F){
    return {F[0], (F[1]-F[0]).det(F[2]-F[0]).unit()};
}
// 点到面的距离，在dir一侧为正，另一侧为负
db disPF(P3 p, P3 p0, P3 dir) { 
	return (p-p0).dot(dir);
}
//点到直线距离
db disLP(P3 k1,P3 k2,P3 q){
    return ((k2-k1).det(q-k1)).abs()/(k2-k1).abs();
}
//直线到直线距离
db disLL(P3 k1,P3 k2,P3 k3,P3 k4){
    P3 dir=(k2-k1).det(k4-k3); if (sgn(dir.abs())==0) return disLP(k1,k2,k3);
    return fabs( dir.unit().dot(k1-k2) );
}
//直线与平面交点
VP getFL(P3 p,P3 dir,P3 k1,P3 k2){
    db a=(k2-p).dot(dir),b=(k1-p).dot(dir),d=a-b;
    if (sgn(fabs(d))==0) return {};
    return {(k1*a-k2*b)/d};
}
//面与面的交线
VP getFF(P3 p1,P3 dir1,P3 p2,P3 dir2){// 返回一条线
    P3 e=dir1.det(dir2),v=dir1.det(e);
    db d=dir2.dot(v); if (sgn(abs(d))==0) return {};
    P3 q=p1+v*dir2.dot(p2-p1)/d; return {q,q+e};
}

//直线与直线交点,平行\重合无输出
VP getLL(P3 pa,P3 pb,P3 qa,P3 qb){
	if( sgn(disLL(pa,pb, qa,qb))!=0 ) return {};
	
	P3 da = (pa-pb).det(qa-qb).unit();
	pair<P3,P3> u = getF({pa ,pa+da, pb});
	return getFL(u.fi,u.se, qa,qb);
}
// 三角形面积
db area(P3 a,P3 b,P3 c){ 
	return ((b-a).det(c-a)).abs();
}
int isMiddle(db a,db m,db b){
	return sgn(a-m)==0||sgn(b-m)==0||(a<m != b<m);
}
// m是否在a构成的b立方体内
int isMiddle(P3 a,P3 m,P3 b){
	return isMiddle(a.x,m.x,b.x)&&isMiddle(a.y,m.y,b.y)&&isMiddle(a.z,m.z,b.z);
}

// 点在线段上
int onSeg(P3 p1,P3 p2,P3 q){
	return sgn((p2-p1).det(q-p1).abs()) == 0 && isMiddle(p1, q, p2);
}
int onSeg_strict(P3 p1,P3 p2,P3 q){
	return sgn((p2-p1).det(q-p1).abs()) == 0 && sgn((q-p1).dot(p1-p2)) * sgn((q-p2).dot(p1-p2)) < 0;
}

// p是否在三角形ps中
int PinT(P3 p,VP ps){
	db a1 = area(p, ps[0], ps[1]);
	db a2 = area(p, ps[1], ps[2]);
	db a3 = area(p, ps[2], ps[0]);
	db a = area(ps[0],ps[1],ps[2]);
	return sgn(a1+a2+a3-a)==0;
}
// 线段AB与三角形ps是否相交, p是交点
int SinT(VP ps,P3 A,P3 B,P3&p){
	P3 v = (ps[1]-ps[0]).det(ps[2]-ps[0]);
	if(sgn(v.dot(B-A))==0) return 0;
	else {
		db t = v.dot(ps[0]-A)/v.dot(B-A);
		if(sgn(t)<0||sgn(t-1)>0)return 0;
		p=A+(B-A)*t;
		return PinT(p,ps);
	}
}
// 三角形与三角形是否相交
int TinT(VP ps,VP qs){
	P3 p;
	rep(i,0,2){
		if(SinT(ps,qs[i],qs[(i+1)%3],p))return 1;
		if(SinT(qs,ps[i],ps[(i+1)%3],p))return 1;
	}
	return 0;
}
// 点p到线段AB距离
db disPS(P3 p, P3 A, P3 B){
	if(A==B)return (p-A).abs();
	P3 v1 = B-A, v2 = p-A, v3 = p-B;
	if(sgn(v1.dot(v2))<0)return v2.abs();
	else if(sgn(v1.dot(v3))>0)return v3.abs();
	else return (v1.det(v2)).abs() / v1.abs();
}
// 线段p到线段q不平行时距离s, p q平行返回0不计算s
int disSS(P3 pa,P3 pb,P3 qa,P3 qb,db&s){
	P3 p1=pa, p2=qa;
	P3 v1 = pb-pa, v2=qb-qa;
	db a1,b1,c1,a2,b2,c2;
	a1=v2.dot(v1), b1=-v1.dot(v1), c1 = -(p2-p1).dot(v1);
	a2=v2.dot(v2), b2=-v1.dot(v2), c2 = -(p2-p1).dot(v2);
	db d = a1*b2 - b1*a2;
	if(sgn(d)==0) return 0;
	db t1,t2;
	t1 = (c1*b2-c2*b1)/d;
	t2 = (a1*c2-a2*c1)/d;
	bool ok1, ok2;
	//q1, q2 两直线最近点对
	P3 q1 = p1 + v1*t2, q2 = p2 + v2*t1;
	if(sgn(disPS(q1, pa,pb))==0)ok1=1;
	else ok1=0;
	if(sgn(disPS(q2, qa,qb))==0)ok2=1;
	else ok2=0;
	if(!ok1&&!ok2)return 0;
	if(ok1 && !ok2) s = disPS(q1,qa,qb);
	else if(ok2&&!ok1) s = disPS(q2,pa,pb);
	else s = (q1-q2).abs();
	return 1;
}
// 点p投影到面(p0, dir)上
P3 getprojF(P3 p,P3 p0,P3 dir){
	return p-dir*((p-p0).dot(dir));
}
// 三角形之间距离
db disTT(VP ps,VP qs){
	db ans = INF;
	if(TinT(ps,qs)) ans=0;
	else {
		rep(i,0,2)rep(j,0,2){
			ans = min(ans, (ps[i]-qs[j]).abs());
		}
		P3 dir1 = ((ps[1]-ps[0]).det(ps[2]-ps[0])).unit();
		P3 tmp;
		rep(i,0,2){
			tmp = getprojF(qs[i], ps[0], dir1);
			if(PinT(tmp, ps)) ans = min(ans, (tmp-qs[i]).abs());
		}
		P3 dir2 = ((qs[1]-qs[0]).det(qs[2]-qs[0])).unit();
		rep(i,0,2){
			tmp = getprojF(ps[i], qs[0], dir2);
			if(PinT(tmp, qs)) ans = min(ans, (tmp-ps[i]).abs());
		}
		rep(i,0,2)rep(j,0,2){
			ans = min(ans, disPS(ps[i], qs[j],qs[(j+1)%3]));
			ans = min(ans, disPS(qs[i], ps[j],ps[(j+1)%3]));
			db s;
			if(disSS(ps[i],ps[(i+1)%3], qs[j],qs[(j+1)%3], s)){
				ans = min(ans, s);
			}
		}
	}
	return ans;
}


// 3D Covex Hull Template
db getV(P3 k1,P3 k2,P3 k3,P3 k4){ // get the Volume
    return (k2-k1).det(k3-k1).dot(k4-k1);
}
db rand_db(){return 1.0*rand()/RAND_MAX;}

VP convexHull2D(VP A,P3 dir){
    P3 x={(db)rand(),(db)rand(),(db)rand()}; x=x.unit();
    x=x.det(dir).unit(); P3 y=x.det(dir).unit();
    P3 vec=dir.unit()*A[0].dot(dir);
    vector<P>B;
    for (int i=0;i<A.size();i++) B.pb(P(A[i].dot(x), A[i].dot(y)));
    B=convexHull(B); A.clear();
    for (int i=0;i<B.size();i++) A.pb(x*B[i].x+y*B[i].y+vec);
    return A;
}

// 使用 ConvexCut3D 求 3D 凸包上三角形面，使用 reduceCH 求出真正的3D凸包
namespace CH3{
    VVP ret; set<pair<int,int> >e;
    int n; VP p,q;
    void wrap(int a,int b){
        if (e.find({a,b})==e.end()){
            int c=-1;
            for (int i=0;i<n;i++) if (i!=a&&i!=b){
                if (c==-1||sgn(getV(q[c],q[a],q[b],q[i]))>0) c=i;
            }
            if (c!=-1){
                ret.pb({p[a],p[b],p[c]});
                e.insert({a,b}); e.insert({b,c}); e.insert({c,a});
                wrap(c,b); wrap(a,c);
            }
        }
    }
    VVP ConvexHull3D(VP _p){
        p=q=_p; n=p.size();
        ret.clear(); e.clear();
        for (auto &i:q) i=i+(P3){rand_db()*1e-4,rand_db()*1e-4,rand_db()*1e-4};
        for (int i=1;i<n;i++) if (q[i].x<q[0].x) swap(p[0],p[i]),swap(q[0],q[i]);
        for (int i=2;i<n;i++) if ((q[i].x-q[0].x)*(q[1].y-q[0].y)>(q[i].y-q[0].y)*(q[1].x-q[0].x)) swap(q[1],q[i]),swap(p[1],p[i]);
        wrap(0,1);
        return ret;
    }
	P3 BaryCenter(VVP A){// 重心
		P3 ans(0,0,0), o(0,0,0);
		db all = 0;
		for(VP nowF: A){
			db vol = getV(o,nowF[0],nowF[1],nowF[2]);
			ans = ans + (o+nowF[0]+nowF[1]+nowF[2])/4.0*vol;
			all += vol;
		}
		return ans/all;
	}
}
VVP reduceCH(VVP A){
    VVP ret; map<P3,VP> M;
    for (VP nowF:A){
        P3 dir=(nowF[1]-nowF[0]).det(nowF[2]-nowF[0]).unit();
        for (P3 k1:nowF) M[dir].pb(k1);
    }
    for (pair<P3,VP> nowF:M) ret.pb(convexHull2D(nowF.se,nowF.fi));
    return ret;
}
// 3D Cut 保留 dot(dir,x-p)>=0 的部分
VVP ConvexCut3D(VVP A,P3 p,P3 dir){
    VVP ret; VP sec;
    for (VP nowF: A){
        int n=nowF.size(); VP ans; int dif=0;
        for (int i=0;i<n;i++){
            int d1=sgn( dir.dot(nowF[i]-p) );
            int d2=sgn( dir.dot(nowF[(i+1)%n]-p) );
            if (d1>=0) ans.pb(nowF[i]);
            if (d1*d2<0){
                P3 q=getFL(p,dir,nowF[i],nowF[(i+1)%n])[0];
                ans.pb(q); sec.pb(q);
            }
            if (d1==0) sec.pb(nowF[i]); else dif=1;
            dif|=(sgn( dir.dot((nowF[(i+1)%n]-nowF[i]).det(nowF[(i+1)%n]-nowF[i])) )==-1);
        }
        if (ans.size()>0&&dif) ret.pb(ans);
    }
    if (sec.size()>0) ret.pb(convexHull2D(sec,dir));
    return ret;
}
// 3D 凸包体积
db vol(VVP A){
    if (A.size()==0) return 0; P3 p=A[0][0]; db ans=0;
    for (VP nowF:A)
        for (int i=2;i<nowF.size();i++)
            ans+=abs(getV(p,nowF[0],nowF[i-1],nowF[i]));
    return ans/6;
}
VVP init(db INF) {
    VVP pss(6,VP(4));
    pss[0][0] = pss[1][0] = pss[2][0] = {-INF, -INF, -INF};
    pss[0][3] = pss[1][1] = pss[5][2] = {-INF, -INF, INF};
    pss[0][1] = pss[2][3] = pss[4][2] = {-INF, INF, -INF};
    pss[0][2] = pss[5][3] = pss[4][1] = {-INF, INF, INF};
    pss[1][3] = pss[2][1] = pss[3][2] = {INF, -INF, -INF};
    pss[1][2] = pss[5][1] = pss[3][3] = {INF, -INF, INF};
    pss[2][2] = pss[4][3] = pss[3][1] = {INF, INF, -INF};
    pss[5][0] = pss[4][0] = pss[3][0] = {INF, INF, INF};
    return pss;
}
// 2D 外接圆圆心
P getOutCir(P p1,P p2,P p3){
	db a1=p2.x-p1.x,b1=p2.y-p1.y,c1=(a1*a1+b1*b1)/2;
	db a2=p3.x-p1.x,b2=p3.y-p1.y,c2=(a2*a2+b2*b2)/2;
    db d=a1*b2-a2*b1;
	return {p1.x+(c1*b2-c2*b1)/d, p1.y+(a1*c2-a2*c1)/d};
}
// 3D 外接圆圆心
P3 getOutCir(P3 a,P3 b,P3 c){
	if( sgn((b-a).det(c-a).abs())==0 ) return {INF,INF,INF};
	VP A = {a,b,c};
    P3 x={(db)rand(),(db)rand(),(db)rand()}; x=x.unit();
	P3 dir = (b-a).det(c-a).unit();
    x=x.det(dir).unit(); P3 y=x.det(dir).unit();
    P3 vec=dir.unit()*A[0].dot(dir);
    vector<P>B;
    for (int i=0;i<A.size();i++) B.pb(P(A[i].dot(x), A[i].dot(y)));
	P C = getOutCir(B[0],B[1],B[2]);
    return x*C.x+y*C.y+vec;
}
// 四面体外接球圆心
P3 getCenter(P3 a,P3 b,P3 c,P3 d){
	P3 d1 = (b-a).det(c-a);
	P3 d2 = (d-a).det(c-a);
	if( sgn(d1.det(d2).abs())==0 && sgn(d1.dot(b-d))==0 ) return P3(INF,INF,INF);
	P3 pa = getOutCir(a,b,c), dira = (b-a).det(c-a).unit();
	P3 pb = getOutCir(a,b,d), dirb = (b-a).det(d-a).unit();
	return getLL(pa, pa+dira, pb, pb+dirb)[0];
}
// 最小包围球
pair<P3,db> min_Ball(VP ps){
	if(ps.size() == 1) return {ps[0],0};
	if(ps.size() == 2) return {(ps[0]+ps[1])/2,(ps[1]-ps[0]).abs()/2};
	P3 tp = getOutCir(ps[0],ps[1],ps[2]);
	if(ps.size() == 3) return {tp,(tp-ps[0]).abs()};
	random_shuffle(ps.begin(),ps.end());
	int n = ps.size();
	P3 o = ps[0]; db r = 0;
	rep(i,1,n-1)if(o.disTo(ps[i])>r+EPS){
		o = ps[i], r = 0;
		rep(j,0,i-1)if(o.disTo(ps[j])>r+EPS){
			o  =(ps[i]+ps[j])/2, r = o.disTo(ps[i]);
			rep(k,0,j-1)if(o.disTo(ps[k])>r+EPS){
				o = getOutCir(ps[i],ps[j],ps[k]), r = o.disTo(ps[i]);
				rep(l,0,k-1)if(o.disTo(ps[l])>r+EPS){
					o = getCenter(ps[i],ps[j],ps[k],ps[l]), r = o.disTo(ps[i]);
				}
			}
		}
	}
	return {o,r};
}

//面切球,球(o,r)在面(p,dir) 的正方向的体积，dir需单位化
db cutFBall(P3 o,P3 r, P3 p,P3 dir){
	P3 u = o + dir * r;
	db h = disPF(u, p,dir);
	if( sgn(h) < 0 ) return 0;
	if(cmp(h,r*2.0)>=0) return 4.0*PI*r*r*r/ 3.0;
	return h*h*(3*r-h)*PI/3.0;
}

int main(){}
